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Abstract 

Brenner has recently proposed modifications to the Navier-Stokes equations that are 
based on theoretical arguments but supported only by experiments having a fairly 
limited range [1,2]. These modifications relate to a diffusion of fluid volume that 
would be significant for flows with high density gradients. So the viscous structure 
of shock waves in gases should provide an excellent test case for this new model. In 
this paper we detail the shock structure problem and propose exponents for the gas 
viscosity-temperature relation based on empirical viscosity data that is independent 
of shock experiments. We then simulate shocks in the range Mach 1.0-12.0 using 
the Navier-Stokes equations, both with and without Brenner's modifications. Initial 
simulations showed Brenner's modifications display unphysical behaviour when the 
coefficient of volume diffusion exceeds the kinematic viscosity. Our subsequent anal- 
yses attribute this behaviour to both an instability to temporal disturbances and a 
spurious phase velocity- frequency relationship. On equating the volume diffusivity 
to the kinematic viscosity, however, we find the results with Brenner's modifica- 
tions are significantly better than those of the standard Navier-Stokes equations, 
and broadly similar to those from the family of extended hydrodynamic models that 
includes the Burnett equations. Brenner's modifications add only two terms to the 
Navier-Stokes equations, and the numerical implementation is much simpler than 
conventional extended hydrodynamic models, particularly in respect of boundary 
conditions. We recommend further investigation and testing on a number of different 
benchmark non-equilibrium flow cases. 
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1 Introduction 



The generally-accepted parameter which indicates the extent to which a local 
region of flowing gas is in thermodynamic equilibrium is the Knudsen number: 

A Ma 

where A is the mean free path of the gas molecules, L is a characteristic length 
of the flow system, the Mach number of the flow Ma = |u|/c with u the flow 
velocity and c the speed of sound, and the Reynolds number Re = p|u|L//i 
with p the mass density and /i the dynamic viscosity. (For high Re flows 
over solid bodies, the characteristic length scale is the boundary layer thick- 
ness, and the denominator on the right hand side of eq. ([T]) is v^Re [3J.) As 
Kn increases, the departure of the gas from local thermodynamic equilib- 
rium increases, and the notion of the g clS clS cl continuum-equilibrium fluid be- 
comes less valid. The range of use of the continuum-equilibrium assumption is 
therefore clearly limited, with the applicability of the classical Navier-Stokes- 
Fourier equations (with standard no-slip boundary conditions) conflned to 
cases where Kn < 0.01, typically. Extended, or modified, hydrodynamics at- 
tempts to extend the range of applicability of the continuum-equilibrium fluid 
model into the so-called 'intermediate-Kn' (or 'transition-continuum') regime 
where 0.01 < Kn < 1. 

Howard Brenner of MIT recently proposed modiflcations to the Navier-Stokes- 
Fourier equations for flows with appreciable density gradients [1]. His theoret- 
ical developments and experimental validations are centred on slow moving 
flows where variations in density are primarily caused by variations in tem- 
perature rather than pressure. Of particular interest to Brenner is the motion 
of particles due to thermal gradients, called 'thermophoresis', which provides 
good, yet fairly limited, supporting evidence for his work. This range of evi- 
dence should be broadened, particularly since his work challenges the funda- 
mentals of conventional fluid dynamics and so demands rigorous validation. 
It is therefore of particular interest to see whether his modiflcations to the 
classical Navier-Stokes-Fourier equations improve their predictive capabilities 
for intermediate-Kn flows. 

In this paper, we investigate the application of Brenner's modifled Navier- 
Stokes equations to the shock structure problem. This is a validation case 

posed extended hydrodynamic models, such as the Burnett equations. Bren- 
ner's modifled equations can be considered an extended hydrodynamic model, 
and their relationship to established models is discussed towards the end of 
this paper. In order to be concise and avoid ambiguity, for the remainder of 
this paper we shall use the expression "Navier-Stokes" to refer to the classi- 
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cal Navier-Stokes- Fourier equation set, and adopt the term "Brenner-Navier- 
Stokes" to refer to Brenner's modified version of these. 



2 The shock structure problem 

The shock structure problem concerns the spatial variation in fluid flow prop- 
erties across a stationary, planar, one-dimensional shock in a monatomic gas. 
We define the flow as moving at a speed u in the positive x-direction, with 
the shock located at x = 0; the upstream conditions at x = — oo are su- 
per/hypersonic and denoted by a subscript '1', downstream conditions at 
X = +00 are denoted by a subscript '2'. While shocks are often modelled 
as discontinuities, their physical properties in fact vary continuously from 
their upstream to their downstream levels over a characteristic distance of 
a few mean free paths because the relaxation times for heat and momentum 
transport are finite. The flow in this shock layer is far from being in local 
thermodynamic equilibrium, typically Kn = 0.2 ~ 0.3 i.e. very much within 
the intermediate-Kn regime. 

Since the Navier-Stokes equations perform poorly at these Kn, the shock 
structure problem is particularly apposite for testing extended hydrodynamic 
models. The problem possesses certain features that make it attractive for 
numerical investigation, particularly if hydrodynamic models with high-order 
derivatives are used: 

(1) there are no solid boundaries to consider; 

(2) the upstream and downstream boundary conditions are clearly defined 
through the Rankine-Hugoniot relations, with all gradients of flow quan- 
tities tending to zero far upstream and downstream of the shock; 

(3) the problem is one- dimensional and steady. 

A monatomic gas possesses no modes of vibration or rotation, cannot dis- 
sociate, and only ionizes at the highest temperatures, so its thermodynamic 
behaviour is generally much simpler than that of polyatomic gases. It is for 
this reason that monatomic gases have generally been preferred as the test 
gas in shock structure experiments and analysis. Resulting data have histori- 
cally been presented in normalised form and the lack of access to the raw data 
presents an opportunity for us to specify the test problem in a nondimension- 
alised form that reduces the pre- and post-processing effortE]: 

• It is convenient to set the upstream temperature Ti = 1 because, for a given 



^ In what follows, the use of a power law viscosity model, described in section [XT} 
is anticipated. 
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yUi, the coefficient of proportionality A is then independent of the exponent 
of s in the power law relation of eq. ([7]). For simplicity we then choose A = 1 
and, hence, fii = 1. 

• It is useful to set Ci = 1 so that the upstream flow speed Ui relates directly 
to the upstream Mach number Mai. The ratio of specific heats, 7, is 5/3 for 
monatomic gases, so the gas constant R = c^/{jT) = 3/5 in our adopted 
nondimensionalised set of units. 

• The final parameters can be set so that a unit length corresponds to the 
upstream Maxwellian mean free path: 

A.n = (2) 

where pi is the upstream pressure. For this test case 16/5^^2717 = 0.99 ~ 1, 
so Aa^i ~ HiCi/pi. By setting pi = 1, the unit length then corresponds 
almost exactly to \mi- 

All normal gradients of p and T are specified as zero at the solution domain 
boundary far downstream. The shock is maintained stationary and fixed within 
the domain by application of the Rankine-Hugoniot velocity relation at the 
downstream boundary, which for our case is 

- = ^ (1 + 3Ma-) . (3) 



The case is initialised with a step change in fields from upstream to down- 
stream at a; = 0. To minimise the time required to reach a steady-state solu- 
tion, the downstream pressure and temperature are initialised using Rankine- 
Hugoniot relations for pressure and temperature, which for our case are: 

^ = i(5Ma?-l) and ^ = Pl^. (4) 
Pi 4 V 1 ; Ti Pi Ml ^ ^ 



The initial and boundary conditions for the actual solution variables, described 
in section m below, are simply derived from those for p, T and u, e.g. initial and 
boundary conditions for p are calculated from p and T using the perfect gas 
law, p = pRT. When we include the second derivative of p through application 
of eqn. f|T^ . below, we specify additionally that the normal gradient of the 
gradient of p is zero. 



3 The viscosity-temperature relation 



One of the main uncertainties in the physical modelling of the shock structure 
problem is the relation between p and T. This is unfortunate because the 
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/i(T) relation has an appreciable effect on the profile of a simulated shock 
— in the extreme case of assuming a constant fi, results of simulations are 
very poor. It could even be argued that the Navier-Stokes equations could be 
made to work for the shock structure problem simply by adjusting the fi(T) 
relation until experimental shock profiles are reproduced. The value of the 
shock structure problem as a good test for hydrodynamic models therefore 
relies on establishing a good /i(T) relation from reliable experimental sources, 
preferably independent of shock data. We therefore review sources to establish 
/i(T) relations for a range of monatomic gases, from which argon is chosen for 
our test problem. 



3. 1 The power law relation 

The viscosity of a perfect rarefied gas is defined through 

= Tip OC TiTlT = -TqUT, (5) 

where n is the gas molecular number density, ti is the collision interval for 
momentum transport, and tq is the collision interval for hard-sphere molecules. 
That the viscosity given in eq. ([5]) is, in fact, purely dependent on temperature 
and not mass density, arises from the nature of the intermolecular force law 
which determines how molecules interact in collision with each other. For 
reasons of simplicity, this force is often modelled, for a given species of molecule 
at a particular temperature, as varying with distance from the molecular centre 
as an inverse power law with coefficient u. In a collision, molecules approach 
each other with a relative speed g and slow to a stop a distance d from each 
other when their kinetic energy is transformed to potential energy in the force 
field, i.e. oc d~'^~^^. With translational temperature a function of the square 
of the molecular velocity, it is then clear that the effective molecular diameter, 
d oc T^^^'^'^^^K The coUisional relaxation time tq is then the mean free path 
(oc divided by the mean molecular velocity (oc T^/^) so that 

nro oc T-^/^c/-^ = T'-\ where s = - + (6) 

Equations (jS]) and yield the well-known relation 

l^ocT' or l^ = AT% (7) 

where A is a constant of proportionality. There are two theoretical limiting 
cases for the intermolecular force law: z/ = oo,s = l/2 corresponds to hard- 
sphere molecules; z/ = 5, s = 1 corresponds to so-called Maxwellian molecules. 
Real molecules generally have a value of u ranging from about 5 to around 15. 
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3.2 Experimental data for monatomic gases 



Increasing the value of the exponent s in eq. ([7]) in a shock structure calculation 
introduces more dissipation, particularly at the high-temperature, downstream 
side of the shock. This additional dissipation acts to smooth out the shock 
layer, increasing its thickness and, in particular, lengthening the downstream 
"tail" in the flow property profiles. 

Therefore, in order to test any hydrodynamic model it is important to set the 
exponent s independently of the shock structure problem under investigation. 
In 1972, Maitland and Smith [21j critically assessed the viscosities of a number 
of gases, obtained from several different sources using a variety of techniques, 
such as the capillary flow method, oscillating disc and rotating cylinder meth- 
ods, and observations of the retardation of an oil drop in free-fall through a 
gas. For the monatomic gases argon, helium and xenon they produced vis- 
cosity data which they estimated to be accurate to 1.5% in the temperature 
range 80-2000K. That upper limit of 2000K is the downstream temperature 
of a shock of Mach 4.3 propagating into a room-temperature gas; hence these 
viscosity data are applicable to shocks of Ma < 4.3. 

In 1958, Amdur and Mason [22] estimated the intermolecular potential at 
higher temperatures from observations of the scattering of high- velocity molec- 
ular beams, and produced tables of the viscosity of gases at temperatures up to 
15000K, corresponding to a shock of Mach 12.5. It is not known how accurate 
these data are; Amdur and Mason estimated that at the higher temperatures 
the error in viscosity could be as much as 10%. 

We have used these data to estimate the value of s (and hence u) in eq. ([6]) for 
different temperature ranges. We have fitted curves of the power law in eq. ([7]) 
by minimising the error in viscosity for the two sets of experimental data. 
More details of this process can be found in [13], but the experimental data 
and best-fit curves for three common monatomic gases are shown in figs. [T] 
and [21 and the corresponding exponents s and u are given in Table [H 

The coefficient u is itself a function of temperature because at higher tempera- 
tures molecules with more energy can penetrate each others' force-fields more 
effectively. Therefore, due to the differences in temperature range and exper- 
imental techniques reported in the Maitland & Smith and Amdur & Mason 
papers, we have considered two ranges of temperature or, equivalently, shock 
Mach number: up to Mach 4.3, and up to Mach 12.5. (It should be noted that 
the power-law fit is not very good for the high temperature data, which is 
itself of unknown accuracy, therefore the fi{T) obtained must be treated with 
caution.) 

In the more limited temperature range of 3500-8500K, Aeschliman and Cam- 
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bel [23] obtained values for argon viscosity to an accuracy of 12% which can 
be represented to within 1% error by a viscosity-temperature exponent of 
s = 0.74. This value compares well with our value of s = 0.76 in Tabled] for 
temperatures in the range 2000-15000K. 

Correlations between direct simulation Monte-Carlo (DSMC) simulations and 
experimental shock density profiles can also provide data for the exponent s, 
particularly at high temperatures [5]IS]l7]l5]l^ll|17j . While it is clearly prefer- 
able in our study of the shock structure problem to use data that are inde- 
pendent of the problem itself, it is worth noting that each of these published 
papers produces a value of u that falls within a range 9 < u < 11, correspond- 
ing to 0.70 < s < 0.75. The value s = 0.72 of Alsmeyer P] is often quoted, 
e.g. recent results of simulations using this value of s agreed with experiment 
within an estimated uncertainty of 5% for temperatures above 2000K [T7] . 

In our present study, three values of the exponent s for argon are therefore 
used: 

• s = 0.68, which is our best-fit for shock Mach numbers up to 4.3; 

• s = 0.76, our best-fit up to Mach 12.5; and 

• s = 0.72, which is the mean of our best-fit values, as well as the commonly 
used value of Alsmeyer [9] that falls in the middle of the range of values 
from DSMC correlations with shock density profiles. As this is the mean 
value, it is used as the 'control' in our study. 

We focus on the power-law form of in this paper in order to discern 

effects on shock structure due to different constitutive models for momentum 
and energy diffusion rather than due to the presumed relationship between 
gas properties. However, we recognise that there are other models available 
for /^(T), e.g. Sutherland's Law, that are generally equivalent to adding a 
weak attractive component to the intermolecular force — which is physically 
more realistic. In our simple power-law model, this attractive force would 
manifest itself as an exponent s that decreases as temperature increases, which 
is generally reported (see, e.g., [21] )■ It is interesting to note, however, that 
our present analysis of experimental data does not appear to bear this out: 
in the case of xenon in Table [H s decreases with increasing temperature, but 
with argon and helium the exponent increases with temperature. 



4 Brenner's modification to the Navier-Stokes equations 

Brenner presents his modification to the Navier-Stokes equations as a change 
in Newton's viscosity law. Before arriving at that discussion, we first express 
the standard governing equations in an Eulerian frame of reference as 
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Conservation of mass 



dp 
dt 



+ V • (u„p) = 
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Conservation of momentum 



dm 



+ V-(u„m) + V-P = 0, 



(9) 



at 



Conservation of total energy 



dE 
'dt 



+ V-(u„E) + V-je + V-(P-u„) =0 



(10) 



where is the mass velocity, the momentum density m = pu^, the total 
energy density E = p{e + |umP/2) with e the specific internal energy, jg 
is the diffusive flux density of internal energy, and P is the diffusive flux 
density of momentum — the familiar stress tensor — defined here as positive 
in compression: P = pi + T, where T is the viscous stress tensor and I the 
unit tensor. 

Based on this sign convention, the constitutive model for a Newtonian fluid 
relates T to the rate of deformation tensor, D, by /i and the bulk viscosity k: 



the symmetric part of a tensor. The deviatoric component of the deformation 



In this constitutive model, eq. (fTTI) . it is generally considered that u is the 
mass velocity that, in the mass continuity equation, relates to a convec- 
tive flux of mass dS • pu^ at an element of surface area dS, or that in the 
Boltzmann equation represents the statistical mean value velocity. However, 
this assumption has recently been questioned by Brenner P|2J who postulated 
that the velocity appearing in Newton's viscosity law should instead be the 
volume velocity u^, so named since it relates to the flux of volume rather than 
mass. 

The distinction between mass and volume flux is perhaps best explained by 
considering a single species fluid at a molecular level. The mass flux through 
dS is the product of the molecular mass and the number of molecules passing 
through dS in one second. There is no net mass flux due to random motions 
of molecules; at a continuum level, there is no diffusive flux of mass, only the 
convective flux already defined. 

To understand volume flux we can consider attributing to each molecule lo- 
cally a microscopic portion of the volume of fluid. The molecular volume is 






for a velocity u, i.e. the overbar indicates 
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transported with the molecule but will change depending on the mass den- 
sity of its surroundings, e.g. the microscopic volume shrinks as the molecule 
moves into a denser region. A convective flux of volume is associated with bulk 
motion of the molecular volumes, and is equivalent to the ratio of convective 
mass flux to mass density, i.e. dS'Um,- If the fluid density varies across dS, 
as a molecule passes through dS there is a change in its associated volume, 
thus a net flux of volume. Random motions can therefore produce a net flux 
of volume, so that there exists a diffusive flux of volume in regions of non-zero 
density gradient. The volume flux dS'u^, therefore represents the total flux 
of volume, comprising the convective flux dS«u,m and a diffusive flux dS«j„, 
where jt, is the diffusive volume flux density, such that 



For a single component fluid undergoing heat transfer, Brenner proposed a 
constitutive equation for j„: 



where ay is termed the 'volume diffusivity' [T]. Exactly how should be quan- 
tified for a given fluid state is an open question. Brenner relates ay directly 
to well-known diffusivity coefficients under some limited conditions; in partic- 
ular, for single component fluids undergoing heat transfer, a^ is the same as 
the thermal diffusivity a = k/pcp, where k is the thermal conductivity and Cp 
is the speciflc heat capacity at constant pressure. 

While the distinction between volume and mass velocities has been made, 
the question remains of why the velocity appearing in Newton's viscosity law 
should be the volume velocity rather than mass velocity. Brenner's justiflca- 
tion is based on limited evidence (e.g. comparison of analytical solutions with 
thermophoresis experiments). A lack of theoretical physical argument could 
therefore lay the hypothesis open to some criticism. However, some support 
for it can be found within the phenomenological GENERIC theory presented 
by Ottinger [25]. First, he demonstrates that the GENERIC formulation ar- 
rives at the standard Navier-Stokes equations when the terms associated with 
mass density in the friction matrix are identically zero. Then, by including 
non-zero terms associated with mass density in the friction matrix, a revised 
set of governing equations is derived that includes two velocities, similar to 
and u„ deflned through eqs. f|T2l) and f|T3l) . What emerges is that the standard 
governing equations have historically ignored mass diffusivity on the basis that 
the diffusive mass flux is zero, while forgetting that there are associated mo- 
mentum and energy fluxes that may not be zero. GENERIC includes these 
momentum and energy fluxes, both of which are entropy producing, making 
the process of mass diffusion irreversible. The ability of mass diffusion to pro- 
duce entropy is, according to Ottinger, something that is missing from the 





3v = 



a^-Vp, 
P 



(13) 
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conventional Navier-Stokes equations. 



Brenner's modification, essentially eq. (fT2l) . can be incorporated into the sys- 
tem of governing fluid equations either by recasting the equations using 
as the convective velocity instead of u^, or by using u^, in the constitutive 
equation for Newton's viscosity law. The former approach has been adopted 
elsewhere [25)126] but here we choose the latter simply so that the Brenner 
modification appears more clearly as a new extended hydrodynamic model, 
rather than a radical change to the governing equations themselves. 

For a monatomic gas, k = 0. Combining eqs. ffTTl) . f|T2|) and f|T3l) yields a 
modified expression for the viscous stress: 



VUm + (VU^ 



2^ 



2/iv|^Vpj 



(14) 



The Brenner approach requires the transport of energy to be similarly modified 
through consideration of the diffusion of internal energy. It is usually assumed 
that the diffusion of internal energy consists solely of heat diffusion, so that 
the diffusive internal energy flux density, j„, is considered synonymous with 
diffusive heat transfer, q = —kVT, according to Fourier's law. However, the 
presence of a diffusive volume flux, of flux density j^, enables energy to be 
transported across a surface by diffusive work transfer of an amount —pjv 
The diffusive internal energy flux density is therefore given by: 

je = q-pj„, (15) 

which, following our argument above, can be re-written in the form: 

je = -kVT - aJVp. (16) 
P 



Equations f|T^ with f|T6|) comprise Brenner's modifications to the classical 
Navier-Stokes- Fourier equations tli2j. We should note that this new fluid 
model has yet to receive either independent theoretical justification or experi- 
mental confirmation. As with any new hypothesis or model it is also subject to 
refinement and re-casting into different forms. However, we use it in this paper 
in the form presented in [T1[2] . without prejudice, to provide an indication of 
its current utility and limitations. 



5 Numerical solution of the governing equations 



Our numerical shock structure solver is developed using the open source Field 
Operation and Manipulation (OpenFOAM) software [27]. Written in C++, 
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OpenFOAM uses finite volume (FV) numerics to solve systems of partial differ- 
ential equations ascribed on any 3-dimensional unstructured mesh of polygonal 
cells. All solvers developed within OpenFOAM are therefore 3-dimensional, 
but can be used for 1- or 2-dimensional problems by the application of par- 
ticular conditions on boundaries lying in the plane of the direction(s) of no 
interest. 

Fluid flow solvers in OpenFOAM are generally developed within an implicit, 
pressure-velocity, iterative solution framework. The solver we developed for 
this work first solves eqs. ([HD, © and flTUl) for p, m and E respectively. The 
equations are treated in a segregated manner: for each equation, terms includ- 
ing the solution variable are, wherever possible, treated implicitly, with other 
terms treated explicitly. All equations include convection of transported vari- 
ables that require a consistent, conservative set of fluxes of u^. After solving 
the sequence of segregated equations for p, m and E, an iterative PISO-style 
method [28] solves an equation for pressure p, derived from the perfect gas 
law, and eqs. ([8]) and iQ, to produce conservative fluxes of m from which the 
fluxes of derived. Finally, m is also corrected from its new fluxes and 

p is corrected from the new solution of p according to the perfect gas law, 
before moving forward to the next time step and returning to the sequence of 
equations for p, m and E. 

Our FV discretisation maintains a compact computational molecule for the 
orthogonal component of the Laplacian terms, which corresponds to Rhie and 
Chow interpolation [2H] in the pressure equation. Both the transported fields 
in convection terms and the fluxes in m are interpolated using limiters from 
the MUSCL total variation diminishing (TVD) scheme [30|l3T] with identical 
limiters used in all convection terms (for p, m and E) to maintain numerical 
consistency. The temporal derivative is discretised using a two-time-level Euler 
implicit scheme. 

We calculated shocks of Mach 1.2, 1.7, 2.2, 2.84, 3.4, 4, 5, 6, 7, 8, 9, 10 and 11 
in order to provide a reasonable distribution of solution points for subsequent 
comparison with results from experiment. Mach 2.84 was specifically chosen 
to coincide with shock profile data communicated to us privately [52]. Simi- 
larly, Mach 8 was chosen to coincide with the published shock profile data of 
Steinhilper |7J, and Mach 9 coincides with a published profile of Alsmeyer [9]. 

A solution domain of 33 Ami was used in all simulations — wide enough to 
contain the entire shock structure comfortably. Our initial results were ob- 
tained using the Navier-Stokes equations with the control viscosity exponent 
s = 0.72. The results for p and T converged on a mesh of 800 cells to within 
1% of the solution extrapolated to an infinitely small mesh size. The results we 
present in this paper were produced with a mesh of 2000 cells, corresponding 
to a mesh size of ~ 0.017Ami- Numerical solutions were executed until they 
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converged to steady-state, at which point the residuals of all equations had 
fallen 5 orders of magnitude from their initial level. 



6 Volume diffusivity 

6.1 Unphysical behaviour when = a 

We followed our initial Navier-Stokes simulations with preliminary simulations 
using the Brenner-Navier-Stokes equations. These simply used a volume diffu- 
sivity = a, which was Brenner's original suggestion (discussed in section H] 
above). However, our simulations do not reach a converged solution with de- 
creasing cell (or mesh) size: at the upstream edge of the temperature profile a 
small undershoot develops at a cell size of O.OGAa/i that increases in magni- 
tude with decreasing cell size. The problem is present in profiles of all solution 
variables but is best illustrated in a plot of Mach number, as shown in fig. [31 At 
best, the level of overshoot at the smallest cell sizes seems unphysical; worse 
is that the overshoot may tend to infinity as the mesh is further refined. 

There is little doubt that the overshoot in Mach number is a consequence 
of Brenner's modification. In subsequent tests we were able to attribute the 
presence of the overshoot to the additional term in the momentum flux, but 
not to the additional term in the energy flux. We therefore postulate that 
the overshoot might be caused by an inappropriately large volume diffusivity, 
particularly since it exceeds the diffusivity associated with the remaining terms 
in the model, the kinematic viscosity u = /i/p, by a factor a* = a^j v = Pr^^ = 
1.5. 

Our search for an alternative value of began by relating the physical pro- 
cess of volume diffusion more closely to mass diffusion, rather than thermal 
diffusion. However, the process of diffusion of mass within a single compo- 
nent fluid, or self- diffusion, occurs at a similar rate to thermal diffusion, with 
theoretical self-diffusivity coefficients Dm of 1.200z/ for hard-sphere molecules 
and 1.534z/ for Maxwellian molecules 03]. The unphysical overshoot in Mach 
number remains for both values of self- diffusion coefficient. The overshoot is 
less pronounced when using the lower value but is still increasing with de- 
screasing cell size even at the smallest cell sizes we tested ( 0.008Aj\/i). Again, 
the amount of overshoot is unphysical and may tend to infinity as the mesh 
is further refined. Only when we reduce the volume diffusivity coefficient to 
ay = V (i.e. when we set a* = 1) does the overshoot disappear. 
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6.2 Investigation of the unphysical behaviour 



It is known that some forms of extended hydrodynamic equations are unsta- 
ble in time to periodic spatial disturbances with wavelengths shorter than a 
critical length that is typically of the order of one mean free path [12] . Such 
instabilities appear in numerical simulations when the mesh is sufficiently fine 
to resolve wavelengths shorter than the critical length, i.e. when the numer- 
ical cell length is below this critical length. The appearance of an overshoot 
below a critical level of mesh refinement in our simulations may indicate a 
similar instability, although the overshoot does appear at a particularly short 
cell length and the solutions do converge to a steady state and so do not 'blow 
up' in time. 

Similarly, some forms of extended hydrodynamic equations, which may be sta- 
ble in time, are actually unstable in space to periodic temporal disturbances 
[M|19|20j . Again, it is unclear that such an instability would produce the 
overshoot behaviour we witnessed in our preliminary calculations. Neverthe- 
less, here we undertake both temporal and spatial stability analyses of the 
Brenner-Navier-Stokes equations in order to investigate the source of unphys- 
ical behaviour. 

Following the procedures described in p^HM] the governing equations from 
section H] are first linearised in 1-dimension to produce the following non- 
dimensionalised perturbation equations: 



dt' ^ 






1 










1 





1 


d(h d 

— H < 

dx' dx' 


a' > 





2 
3 







¥ 



(17) 



where t' and x' are nondimensionalised time and distance, respectively, = 
{p' m' T'Y^ is the vector of nondimensionalised density, flow speed and tem- 
perature, and a' and g' are nondimensionalised momentum and heat fluxes 
respectively. From eqs. (1141) and (I16p . these momentum and energy fluxes are, 
respectively, 

4 9m' 4 

a = — fy* 

3 ox' 

and 

15 dT' 



3 dx'"^ 



a 



4 dx' " dx' 
of the form 



(19) 



We assume a solution to eq. ([T] 

</) = 0exp{i(cut'-A;x')}, (20) 
where is the amplitude of the wave, uj is its frequency and k its propagation 
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constant. Equations (fT71) to (120!) can be combined to produce a set of linear 
algebraic equations of the form 



A{u,k)(l) = 0, 



(21) 



for which non-trivial solutions require 



det[A{uj,k)] = 0. 



(22) 



For the Brenner-Navier- Stokes equations, eq. (122]) yields the following charac- 
teristic equation: 



6icj3 + 23A;^cu^ - [lOA;^ + (20 + 8a*)k^]iuj - [(15 - Aa*)k^ + 20a*k^] = 0. (23) 



If a disturbance in space is considered as an initial-value problem, k is real 
and u = ujr + iuji is complex. The form of eq. (!20|) indicates that stability 
then requires cuj > 0. If a disturbance in time is considered as a problem of 
signalling from the boundary, u is real and k = kr + iki is complex. For a wave 
travelling in the positive x direction, kr > 0, and stability then requires that 
ki < 0. For a wave travelling in the negative x direction, the converse is true: 
kr < and stability requires fc, > 0. 

We examine temporal stability by solving eq. (1251) numerically for uj for values 
of k in the range < k < oo. Trajectories of u are plotted in the complex plane 
in fig. m Two sets of trajectories are plotted: those for = a (corresponding 
to a* = Pr~^ = 1.5) and those for = i/ (for which a* = 1.0). In both 
cases the trajectories all lie within the region Ui > 0, indicating stability for 
all k. This confirms, as expected, that the observed overshoot is not caused 
by temporal instability. 

We then turn to examine spatial stability by solving eq. (12^ numerically for 
k for values of uj in the range < u; < oo. Trajectories of k are plotted in 
the complex plane in fig. [5] for both = u and = a. When a^^ = u the 
trajectories do not violate the stability condition. However, when = a, 
the inset graph shows one trajectory start from = 0, fcj. = at fcj ~ 0.55, 
enter the unstable region {kr > 0,ki > 0}, and exit into the stable region 
{kr < 0,ki > 0} by crossing the kr = axis at ki ^ 0.56. Thus, the stability 
condition is clearly violated for small u. 

Subsequently we examined trajectories for a number of different a* and we 
found that the equations are unstable for a* > 1.45, which suggests a potential 
problem for some intuitive choices of a^, such as a and Dm- However, this 
result does not really explain the cause of the overshoot in Mach number, 
since the overshoot only disappears when a* falls to unity, i.e. considerably 
lower than 1.45. 
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Unphysical behaviour can also be observed by examining the phase velocity: 

Figure [6] shows the phase velocity, normalised by the speed of sound in the 
nondimensionalised form in the perturbation equations (cq = a/7), for a* = 1.0 
and a* = 1.2. For both a*, results for mode 1 are superimposed and correspond 
to the propagation of sound. The mode 2 results correspond to the diffusive 
transport of heat and results for both a* are very similar. However, there is 
a marked difference between results for the two cases for mode 3, relating to 
higher-order diffusive transport. The a* = 1.0 results are similar to those of 
other extended hydrodynamic models, e.g. the super-Burnett equations [34j, 
beginning at a moderate speed, Vph/co = 2.34 at = 0, before increasing 
steadily with increasing uj. The a* = 1.2 results are, however, unusual: the 
phase velocity at a; = 0, i.e. Vph/co = 3.99, is high in comparison to other 
hydrodynamic models [34J. The phase velocity also decreases initially with 
increasing uj, before passing through a minimum and increasing thereafter. 
The high initial phase velocity seems anomalous, and rises to extraordinary 
levels for higher a*, e.g. if a* = 1.5, Vph/co = 190.1 at a; = 0. The initial 
decrease in phase velocity with increasing u may allow low frequency waves 
upstream of the shock to overtake slower, higher frequency waves within the 
shock, creating counter-dispersion at the upstream end of the shock. The initial 
decrease in phase velocity disappears only when a* falls below ~ 1.11, a level 
quite close to that at which we find the unphysical overshoot disappears. 

To summarise, our results show unphysical behaviour for the Brenner-Navier- 
Stokes equations when a* = 1.5. A spatial stability test confirms the equations 
are unstable to temporal disturbances when a* > 1.45. Plots of phase velocity 
raise further questions about the physical nature of the solutions when a* > 
1.11. Taken together, this casts doubt both on Brenner's proposed = a and 
on the apparently natural choice of = Dm- The overshoot in Mach number 
disappears when a* = 1.0, i.e. = u. We therefore adopt this model for 
in the Brenner-Navier-Stokes equations for the remainder of this paper. 



7 Results and comparison with experiment 



7. 1 Shock profiles 



We prefer, where possible, to compare results with actual experiments rather 
than DSMC simulations, since the latter requires certain assumptions relating 
to the form of the intermolecular force law. We therefore first make comparison 
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with actual measured data for the variation of density through the shock layer 



Figure [7] shows the normalised variation of density and temperature, p* and 
T* respectively, through an argon gas shock of Mach 2.84 calculated using 
the Navier-Stokes and Brenner-Navier- Stokes equations with s = 0.72. The 
experimental density profile [32] is also shown. It is clear that the shock layer 
predicted by the conventional Navier-Stokes equations is too thin, whereas 
the Brenner-Navier-Stokes equations produce good agreement with the exper- 
imental data. The main region of disparity is upstream of the shock layer (left 
hand side in the figure) where the experimental data trails out and is flatter 
than the prediction. 

Similarly, fig. [8] shows the predicted profiles for a Mach 8.0 shock compared 
with experimental density data [7|. Again, the Navier-Stokes equations pro- 
duce a shock profile which is too thin when compared with experiment. How- 
ever, the Brenner-Navier-Stokes equations produce excellent agreement in the 
central and downstream shock regions (p* > 0.2), only in the upstream region 
is the predicted profile sharper than the experimental data shows — just as 
in the Mach 2.84 case. 

Figure [H] shows results for a Mach 9.0 shock In this case our observations 
are very similar to those we make about the Mach 8.0 shock profile, above; 
this figure is included here for completeness. 

7.2 Inverse density thickness 

Apart from direct comparison of calculated and experimental shock profiles, 
there are other shock parameters for which experimental and/or independent 
numerical data is available. The principal parameter is the non-dimensional 
shock inverse density thickness, defined as: 

1 ^ ^M^|y^|^^^_ ^25) 

P2- Pl 

In the absence of an a priori characteristic length scale in an unconfined flow, 
the definition of Kn requires a characteristic dimension of a flow structure, in 
this case the actual thickness of the shock layer itself. Therefore has the 
interesting feature that it represents Kn for the shock structure casa I. 

Alsmeyer [9] reported the most comprehensive collection of experimental shock 

^ While this identification then indicates, as we see below, that shocks generally 
have such a high overall Kn that any hydrodynamic model should fail, we can still 
assess extended hydrodynamic equations for their usefulness as engineering models. 
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data, comprising previously-published work as well as his own new results. 
Figure [TO] shows for argon shocks up to Mach 11, with experimental data 
collated from a number of sources [7ll9ii32j . The Navier-Stokes equations, with 
s = 0.72, predict shocks of approximately half the measured thicknesses over 
the entire Mach number range. As indicates Kn, this poor agreement is 
expected because we can see that Kn ~ 0.2 - 0.25 over most of this Mach 
number range, so the Navier-Stokes equations are beyond their effective range 
of application. However, the results from the Brenner-Navier-Stokes equations 
closely match experiment, with moderate sensitivity to the choice of viscosity- 
temperature exponent: using 0.72 (the control value) and 0.76 (our best-fit 
value for temperatures equivalent to shocks up to Mach 12.5), the results fall 
within the limits of experimental scatter; using s = 0.68, the results stray 
slightly outside the scatter of experimental data just before they reach the 
exponent's limit of applicability at Mach 4.3. With the results for 0.72 at the 
higher end of the experimental scatter, and results for 0.76 at the lower end, 
we estimate that an exponent of s ~ 0.74 would produce the best agreement 
with the experimental results. 



7.3 Density asymmetry quotient 



Agreement of predicted and experimental shock inverse density thicknesses is 
not the only measure of the success of a new model. As L^^ depends on the 
density gradient at the profile midpoint alone, it does not express anything 
about the overall shape of the profile. If is used as the sole parameter to 
describe the shock it could be concluded that the Brenner-Navier-Stokes equa- 
tions tested here, and many other models previously published, have excellent 
predictive capability. However, the shock profiles in figs. [7] and [8] show there 
are differences between simulation and experiment, in particular relating to 
the flatter region upstream of the profile that is observed experimentally. 

Therefore, a second parameter which should be used to describe the shock 
profile, and for which experimental data is available, is the density asymmetry 
quotient Qp. This is a measure of how skewed the shock density profile is 
relative to its midpoint. It is defined for a 1-dimensional profile of normalised 
density, p*, centred at = 0.5 on x = 0, as 

A symmetric shock would consequently have Qp = 1, but real shock waves 
are not completely symmetrical about their midpoint. First, their 'bulk' form 
is generally skewed somewhat towards the downstream. Then, the aforemen- 
tioned flattened, diffusive region, that extends upstream of the shock profile, 
tends to increase Qp. Figure [TT] shows experimental data compiled by Alsmeyer 
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[9] in which Qp increases fairly hnearly from ~ 0.9 at around Mach 1.5, through 
unity at around Mach 2.3, to ~ 1.15 at Mach 9. The bulk form therefore cor- 
responds to Qp ~ 0.9 and the upstream flattened region accounts for a further 
increase in Qp, of up to 0.25 at Mach 9. 

Results from the Navier-Stokes equations do not agree well with experimental 
data: the bulk form is skewed towards the upstream side so that Qp > 1 
even at the lowest Mach numbers and the skewness further increases with 
Mach number (apparently by sharpening of the profile downstream rather 
than flattening upstream) so that by Mach 4, Qp f» 1.4, compared to ~ 1.03 
from experiment. 

We find the Brenner-Navier-Stokes equations predict the bulk form of the 
profile very well, predicting Qp 0.9 at low Mach numbers. As discussed 
in section I7.H it does not capture the flattened region upstream and so the 
departure from experimental data increases with Mach number. 



7.4 Temperature- density separation 



The final shock structure parameter is the temperature-density separation, 
5tpi measured between the midpoints of the respective normalised profiles. In 
a shock, the density rises from its upstream value to its downstream value 
behind the temperature, due to the finite relaxation times for momentum and 
energy transport; a good hydrodynamic model should resolve this spatial lag 
accurately. However, experimental data for this parameter is scarce due to the 
difficulty in measuring temperature profiles, so independent DSMC data [TT] 
is usually taken for comparison. 

Figure [12] shows the temperature and density profiles for a Mach 11 shock 
calculated using DSMC [TT], and the Navier-Stokes and the Brenner-Navier- 
Stokes models. As in our earlier comparisons in section 17.1^ the Brenner- 
Navier-Stokes equations produce profiles that are much sharper in the up- 
stream region of the shock than those from DSMC results. 

Figure [12] compares DSMC data for dxp over a range of shock Mach numbers 
with results from our simulations. The DSMC data show an increase in 5tp 
from ~ 1.5Aa/i at Mach 1.5 to ~ 2.9AAfi at Mach 8. The Navier-Stokes and 
Brenner-Navier-Stokes equations increasingly under-predict 6tp with increas- 
ing Mach number, although the Brenner-Navier-Stokes equations generally 
perform a little better over the Mach number range. 
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7.5 Very strong shocks 



The inverse density thickness of extremely strong shocks is a useful additional 
comparison for any proposed hydrodynamic model. Narasimha and Das [10] 
examined the solution of the Boltzmann equation for an infinitely strong shock 
(a more recent treatment is in [U]), modelling the upstream flow as a molec- 
ular beam with a distribution function of the form f{x = — oo) = ni5(ui), 
where 6 is the Kronecker delta function. The shock layer may then be treated 
as a device for converting this beam function into a downstream Maxwellian 
distribution function. The distribution function in the shock layer can be ex- 
pressed as a linear combination of the two extremal distribution functions, a 
method similar to the bimodal method of |4j . 

Using an expansion parameter that measures the departure of the distribution 
function in the shock wave from that outlined in the previous paragraph, an 
infinite series of ordinary differential equations is obtained for the shock thick- 
ness. This series rapidly converges, and a solution of the first seven equations 
of the set yields a predicted shock thickness of 6.7 Xm2 (which is written in 
terms of the downstream Maxwellian mean free path version of eq. [2]) . When 
this is converted into the of eq. (!25l) . the inverse density thickness for a 
shock with a downstream Xm2 equivalent to that for a shock of Mach 100 is 
predicted to be 0.076. 

Our calculations for shocks of Mach 100 give = 0.156 for the Navier- 
Stokes equations with s = 0.72. The results with the Brenner-Navier-Stokes 
equations are = 0.091 for s = 0.72, and = 0.066 for s = 0.76. These 
values for straddle the solution from the molecular beam analysis. Further 
simulations with successive adjustments to s gave a precise match in for 
s = 0.742, which is in agreement with the exponent estimated in section 
to produce the best agreement with experiment over the range Mach 1-11. 



8 Discussion and conclusions 

The Navier-Stokes equations are robust and accurate over a wide range of Kn 
— surprisingly so, given some of the relatively narrow axioms on which they 
depend (i.e. the continuum- fluid and local-equilibrium requirements). Such a 
good fluids engineering model is difficult to relinquish, even when flow systems 
well beyond its range of applicability are considered [35]. However, it is clear 
from the results in section [7] that the Navier-Stokes equations fail in nearly 
every respect in predicting correct shock structures above about Mach 2 (or, 
equivalently, for intermediate-Kn flows). 
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While it is important not to draw strong conclusions based on just one test 
case, our results are generally encouraging for the Brenner-Navier-Stokes equa- 
tions. This modified model is significantly better at reproducing the trends in 
the experimental and DSMC data, and in the case of the inverse density thick- 
ness delivers an excellent match. It is only the more detailed features of the 
shock profile that Brenner's model seems unable to reproduce. 

First, it does not predict the fiattened upstream region, as discussed in sec- 
tion l7.3[ In this regard, a major advantage of DSMC as a technique for simulat- 
ing intermediate-Kn fiows in general is its ability to produce non-Maxwellian 
velocity distributions, that may also differ in directions parallel and perpen- 
dicular to the flow. It is not clear that hydrodynamic models will be able 
to properly incorporate this physics, and certainly the problems the Brenner- 
Navier-Stokes equations have in capturing the upstream shock region properly 
is related to the fact that in this region the velocity distribution function is a 
non-Maxwellian combination of fast, cold upstream molecules and slower, hot 
molecules that have diffused from downstream regions. 

The second feature is also related to this distribution function problem: bi- 
modal methods (see, e.g., |15j) for a hard-sphere gas predict a downstream 
temperature overshoot of around 1%, which is confirmed by careful DSMC 
simulations. There are no downstream overshoots predicted in any of the 
Brenner-Navier-Stokes shock simulations. 

While some of these features can be obtained using certain extended hydro- 
dynamic models that are formally O(Kn^) (see, e.g., pT|IT3] ). this is at a 
cost: there are known problems of physical stability, and the numerical im- 
plementation is difficult due to the large number of additional non-linear and 
high order derivatives. The Brenner modification does not suffer so much from 
these problems, having only a single additional term in each of the momentum 
and energy conservation equations. That the adoption of these terms provides 
a substantial improvement in predicted results raises the question of whether 
this model can compensate, in part at least, for increased non-local-equilibrium 
in the gas, or whether this agreement is coincidental. 

Brenner proposed his modifications partly to understand how some effects that 
are traditionally thought of as becoming important only in a flow approach- 
ing the intermediate-Kn range, e.g. slip at solid bounding surfaces, can be 
encompassed in a model which still retains its essential 0(Kn) character [T||2] . 
He shows that the form (if not the exact coefficients, except for a particular 
molecular model) of two particular terms that appear in the Burnett consti- 
tutive model for T are in fact encompassed by the additional term in p in his 
modified Newtonian T of eqn. f|T^ . While all the stress terms in the Burnett 
equations are formally O(Kn^), under some circumstances these two terms can 
be of similar magnitude to those of 0(Kn) i.e. the same order of accuracy as 
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the Navier-Stokes equations. If the issue of the correct model for the volume 
diffusivity, a^, can be resolved then the Brenner-Navier-Stokes equations may 
provide a simple alternative to the family of extended hydrodynamic models 
that includes those of Burnett, Grad, etc., producing reasonably accurate so- 
lutions of intermediate-Kn flows at a modest computational cost. While it is 
known that the classical Burnett equations do not satisfy the second law of 
thermodynamics, truncated or extended forms of the equations can be con- 
structed that do pnil2|13j . The fact that the Brenner-Navier-Stokes equations 
are less prone to both numerical instability and unphysical solution may indi- 
cate that thermodynamic consistency is less of a problem with these equations 
than with more complex extended hydrodynamics models. 

We recognise that it is not reasonable to rely on one benchmark case to de- 
cide the value of the Brenner-Navier-Stokes equations, or any other extended 
hydrodynamic model (or its associated boundary conditions). Equation ([T]) 
indicates there are three distinct categories of near-equilibrium flows: 

A. Ma = 0(1), Re oo, typical of super- and hypersonic flows; 

B. Ma — » 0, Re = 0(1), typical of flows in micro- or nano-systems; 

C. Ma — > 0, Re — > oo, typical of incompressible turbulent boundary layer 



As Kn vanishes more quickly for flows in category C than in categories A and 
B, departures from local equilibrium in category C flows are not as signiflcant 
as in categories A and B flows. This paper has addressed a category A flow 
in which the boundary conditions are not in doubt, but benchmark cases 
for models of intermediate-Kn flows generally require additional boundary 
conditions, usually at solid surfaces, the speciflcation of which is one of the 
outstanding problems in hydrodynamic approaches to rarefled gas dynamics. 
Setting aside the boundary condition problem, however, we can propose a 
number of benchmark cases in categories A and B that any new hydrodynamic 
model for rarefled flows should be tested against: 

• the shock structure problem, as outlined in this paper (including compar- 
isons of Qp, and the thickness of Ma = oo shocks, in addition to the 



• the nonlinearity of the thermal and momentum Knudsen layers (the region 
0(A) close to sohd surfaces); 

• the 'Knudsen paradox' in Poiseuille flow, i.e. the minimum in the mass flow 
rate at around Kn ^ 1, as well as bimodality in the temperature proflle; 

• drag coefficients and heat transfer in: flow around a sphere, flow around a 
cylinder and Couette flow; 

• variation of skin friction on cones in supersonic flow; 

• base pressures on cone-cylinder conflgurations in supersonic flow as the 
Knudsen layer extends far into the wake region; 



flows. 
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• thermophoresis, i.e. the force on particles suspended in a rarefied gas be- 
tween two parallel plates of different temperature; 

• dispersion and absorption of ultrasonic sound waves. 

Tfiis list is neither exclusive nor comprehensive; we are sure that other good 

benchmark cases could be added to it. The caveat is that in most cases ex- 
perimental data is extremely sparse and unreliable, and unfortunately much 
reliance still needs to be placed on comparison with independent DSMC or 
other molecular dynamics simulations as 'experimental analogues'. 
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Tables 





Argon 


Helium 


Xenon 


Applicable range 


s 




s 


u 


s 




up to Mach 4.3 


0.68 


12.0 


0.71 


10.7 


0.77 


8.5 


Mach 4.4-12.5 


0.76 


8.8 


0.83 


7.2 


0.72 


10.3 



Table 1 

Collated experimental values of s and u for argon, helium and xenon gases. 
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Figures 




10 ^ < < < < ^ 

400 800 1200 1600 2000 



Temperature, T [K] 

Fig. 1. Experimental viscosity versus temperature data for intermediate tempera- 
tures for argon, helium and xenon, fitted to a power- law: oc T*. 
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Fig. 2. Experimental viscosity versus temperature data for high temperatures for 
argon, hehum and xenon, fitted to a power-law: /x oc T*. 
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Fig. 3. Mach number profiles for decreasing cell sizes, assuming a„ = a. 
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Fig. 4. Temporal stability of Brenner-Navier-Stokes equations (arrows indicate di- 
rection of increasing k; grey shaded area indicates region of instability). 
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Fig. 5. Spatial stability of Brcnncr-Navier-Stokes equations (arrows indicate direc- 
tion of increasing w; grey shaded areas indicate regions of instability). 
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Fig. 6. Normalised phase velocity for the Brenner-Navier-Stokes equations. 




.6 -5 -4 -3 -2 -1 1 2 3 4 5 6 



Distance through shock [Ami] 

7. Simulated and experimental profiles of a Mach 2.84 stationary shock; 
0.72. 
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Fig. 8. Simulated and experimental profiles of a Mach 8.0 stationary shock; s = 0.72. 
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Fig. 9. Simulated and experimental profiles of a Mach 9.0 stationary shock; s = 0.72. 
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Fig. 10. Simulated and experimental inverse density thickness (Lp^) data, versus 
shock Mach number; for various values of the exponent s. 
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Fig. 11. Simulated and experimental asymmetry quotient (Qp) data, versus shock 
Mach number; s = 0.72. 
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Fig. 12. Simulated and DSMC profiles of a Mach 11 stationary shock; s = 
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Fig. 13. Simulated and independent DSMC temperature-density separation 
data, versus Mach number; s = 0.72. 
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